home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / fx_root.pro < prev    next >
Text File  |  1997-07-08  |  6KB  |  153 lines
  1. ;$Id: fx_root.pro,v 1.4 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       FX_ROOT
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes real and complex roots (zeros) of
  11. ;       a univariate nonlinear function.
  12. ;
  13. ; CATEGORY:
  14. ;       Nonlinear Equations/Root Finding
  15. ;
  16. ; CALLING SEQUENCE:
  17. ;       Result = FX_ROOT(X, Func)
  18. ;
  19. ; INPUTS:
  20. ;       X :      A 3-element initial guess vector of type real or complex.
  21. ;                Real initial guesses may result in real or complex roots.
  22. ;                Complex initial guesses will result in complex roots.
  23. ;
  24. ;       Func:    A scalar string specifying the name of a user-supplied IDL
  25. ;                function that defines the univariate nonlinear function.
  26. ;                This function must accept the vector argument X.
  27. ;
  28. ; KEYWORD PARAMETERS:
  29. ;       DOUBLE:  If set to a non-zero value, computations are done in
  30. ;                double precision arithmetic. 
  31. ;
  32. ;       ITMAX:   Set this keyword to specify the maximum number of iterations
  33. ;                The default is 100.
  34. ;       
  35. ;       STOP:    Set this keyword to specify the stopping criterion used to
  36. ;                judge the accuracy of a computed root, r(k). 
  37. ;                STOP = 0 implements an absolute error criterion between two
  38. ;                successively-computed roots, |r(k) - r(k+1)|.
  39. ;                STOP = 1 implements a functional error criterion at the 
  40. ;                current root, |Func(r(k))|. The default is 0.
  41. ;
  42. ;       TOL:     Set this keyword to specify the stopping error tolerance.
  43. ;                If the STOP keyword is set to 0, the algorithm stops when
  44. ;                |x(k) - x(k+1)| < TOL.
  45. ;                If the STOP keyword is set to 1, the algorithm stops when 
  46. ;                |Func(x(k))| < TOL. The default is 1.0e-4.
  47. ;
  48. ; EXAMPLE:
  49. ;       Define an IDL function named FUNC.
  50. ;         function FUNC, x
  51. ;           return, exp(sin(x)^2 + cos(x)^2 - 1) - 1 
  52. ;         end 
  53. ;
  54. ;       Define a real 3-element initial guess vector.
  55. ;         x = [0.0, -!pi/2, !pi]
  56. ;
  57. ;       Compute a root of the function using double-precision arithmetic.
  58. ;         root = FX_ROOT(x, 'FUNC', /double)
  59. ;
  60. ;       Check the accuracy of the computed root.
  61. ;         print, exp(sin(root)^2 + cos(root)^2 - 1) - 1
  62. ;
  63. ;       Define a complex 3-element initial guess vector.
  64. ;         x = [complex(-!pi/3, 0), complex(0, !pi), complex(0, -!pi/6)]
  65. ;
  66. ;       Compute a root of the function.
  67. ;         root = FX_ROOT(x, 'FUNC')
  68. ;
  69. ;       Check the accuracy of the computed root.
  70. ;         print, exp(sin(root)^2 + cos(root)^2 - 1) - 1
  71. ;
  72. ; PROCEDURE:
  73. ;       FX_ROOT implements an optimal Muller's method using complex 
  74. ;       arithmetic only when necessary.
  75. ;
  76. ; REFERENCE:
  77. ;       Numerical Recipes, The Art of Scientific Computing (Second Edition)
  78. ;       Cambridge University Press
  79. ;       ISBN 0-521-43108-5
  80. ;
  81. ; MODIFICATION HISTORY:
  82. ;       Written by:  GGS, RSI, March 1994
  83. ;       Modified:    GGS, RSI, September 1994
  84. ;                    Added support for double-precision complex inputs.
  85. ;-
  86.  
  87. function fx_root, xi, func, double = double, itmax = itmax, $
  88.                             stop = stop, tol = tol  
  89.  
  90.   on_error, 2 ;Return to caller if error occurs.
  91.  
  92.   x = xi + 0.0 ;Create an internal floating-point variable, x.
  93.   sx = size(x)
  94.   if sx[1] ne 3 then $
  95.     message, 'x must be a 3-element initial guess vector.'
  96.  
  97.   ;Initialize keyword parameters.
  98.   if keyword_set(double) ne 0 then begin 
  99.     if sx[2] eq 4 or sx[2] eq 5 then x = x + 0.0d $
  100.     else x = dcomplex(x)
  101.   endif
  102.   if keyword_set(itmax)  eq 0 then itmax = 100
  103.   if keyword_set(stop)   eq 0 then stop = 0
  104.   if keyword_set(tol)    eq 0 then tol = 1.0e-4
  105.  
  106.   ;Initialize stopping criterion and iteration count.
  107.   cond = 0  &  it = 0
  108.  
  109.   ;Begin to iteratively compute a root of the nonlinear function.
  110.   while (it lt itmax and cond ne 1) do begin
  111.     q = (x[2] - x[1])/(x[1] - x[0])
  112.     pls = (1 + q)
  113.     f = call_function(func, x)
  114.     a = q*f[2] - q*pls*f[1] + q^2*f[0]
  115.     b = (2*q+1)*f[2] - pls^2*f[1] + q^2*f[0]
  116.     c = pls*f[2]
  117.     disc = b^2 - 4*a*c
  118.     roc = size(disc)  ;Real or complex discriminant?
  119.     if roc[1] ne 6 and roc[1] ne 9 then begin  ;Proceed toward real root.   
  120.       if disc lt 0 then begin  ;Switch to complex root.
  121.         ;Single-precision complex.
  122.         if keyword_set(double) eq 0 and sx[2] ne 9 then begin
  123.           r0 = b + complex(0, sqrt(abs(disc)))
  124.           r1 = b - complex(0, sqrt(abs(disc)))
  125.         endif else begin ;Double-precision complex.
  126.           r0 = b + dcomplex(0, sqrt(abs(disc)))
  127.           r1 = b - dcomplex(0, sqrt(abs(disc)))
  128.         endelse
  129.         if abs(r0) gt abs(r1) then div = r0 $  ;Maximum modulus.
  130.           else div = r1
  131.        endif else $
  132.         div = max([abs(b + sqrt(disc)), abs(b - sqrt(disc))]) ;Real root. 
  133.     endif else begin  ;Proceed toward complex root.
  134.       c0 = b + sqrt(disc)
  135.       c1 = b - sqrt(disc)
  136.       if abs(c0) gt abs(c1) then div = c0 $  ;Maximum modulus.
  137.         else div = c1
  138.     endelse
  139.     root = x[2] - (x[2] - x[1]) * (2 * c/div) 
  140.     ;Absolute error tolerance.
  141.     if stop eq 0 and abs(root - x[2]) le tol then cond = 1 $ 
  142.     else $  
  143.     ;Functional error tolerance. 
  144.     if stop ne 0 and abs(call_function(func, root)) le tol then cond = 1
  145.     x = [x[1], x[2], root] 
  146.     it = it + 1
  147.   endwhile
  148.   if it ge itmax and cond eq 0 then $
  149.     message, 'Algorithm failed to converge within given parameters.'
  150.   return, root
  151.  
  152. end
  153.